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ABSTRACT 


This report describes the research performed during the period July 
1989 - July 1990 under a joint research program between the 
Mechanical Engineering Department of the University of Maryland and 
the Center for Fire Research of the National Institute of Standards 
and Technology. The research is conducted by Graduate Research 
Assistants of the ME Department under the joint supervision of Dr. 


ai Marzo (UMCP) and of Dr. Evans (CFR - NIST). 


A new experimental set-up for the study of dropwise evaporation in 
a radiant heat transfer field has been designed, constructed and 
tested. The various issues of concern such as: steady state solid 
temperature distribution, radiant heater design and configuration, 
infrared background noise and post test data manipulation are 


outlined. 


The formulation of a model for the prediction of the cooling 
induced by an evaporating droplet impinging a semi-infinite solid 
is the subject of this report. The thermal interactions during the 
evaporation of a liquid droplet deposited on a low conductivity 
semi-infinite solid are complex because the evaporative process is 
coupled to the solid intense local cooling. A combined Boundary 
Element Method (BEM) and Control Volume Method (CVM) has been used 
to solve this complex numerical problem. The results for both high 


and low thermal conductivity materials are in excellent agreement 
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with the experimental findings. Detailed comparison of the surface 
temperature distributions on solid Macor detected with infrared 
thermography are also performed to demonstrate the accuracy of the 


computational method. 
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1. INTRODUCTION 


1.1 Background 


The long term objective of the study of droplet-solid interaction 
is to obtain information applicable to the extinguishment of fire 
through a droplet array (e.g. spray). The solids of concern 
include low thermal conductivity materials, typical of fire 


applications. 


This report describes the research performed during the period July 
1989 - July 1990 under the joint research program conducted by the 
Center for Fire Research (NIST) and the Mechanical Engineering 
Department (University of Maryland) addressing the Transient 
Cooling of a Hot Surface by Droplet Evaporation. iret yO. Tio 


research program was initiated in January 1985. 


The research described hereafter constitutes portion of an 
extensive research program aimed at developing accurate droplet 
cooling models of burning solid fuel surfaces. Many studies have 
been performed to quantify the vaporization process of both single 
and multi droplet arrays impacting on hot surfaces. For the 
studies found in the published literature, the full span of the 
droplet vaporization processes are usually reported. These would 
include evaporation, nucleate boiling, film boiling and Leidenfrost 


transition. This present research is more limited in the span of 


vaporization process studied, being only concerned with the 
evaporation of single droplets on a hot surface. Limited 


experiments are also performed in the nucleate boiling range. 


Several important results were obtained in the first years of 
research. In particular, the modelling of the boundary condition 
at the liquid-vapor interface (at the droplet exposed surface) was 
validated with the data collected for water droplets evaporating on 
an aluminum block (diMarzo 1986a, 1986b, 1988). The simple model 
that describes the cooling effect induced in the aluminum by the 
evaporating droplet (diMarzo 1986a, 1987b, 1987c) provided the 
input for the formulation of a preliminary multi-droplet model 


(diMarzo 1987a). 


An extensive experimental data base was also collected for water 
droplet evaporating on Macor (a glass like material able to 
withstand strong thermal stresses). An infrared thermographic 
technique has been developed to monitor the temperature 
distribution on the surface of Macor during the evaporative 
process. These data are used to validate the more complex model 
required to describe the coupled liquid-solid interaction for the 


Macor case. 


1.2 Overview 


The rationale for the more complex coupled model has been described 


at length in the previous reports. The low thermal conductivity of 
Macor invalidates the assumption of constant uniform temperature 
under the liquid droplet which was crucial in developing the simple 
model for the case of aluminun. This fact requires that the 
droplet evaporative process be solved jointly with the cooling 


transient in the solid. 


A combined Boundary Element Method (BEM) was proposed by Dr. Baum 
(CFR - NIST). This method is used to integrate the solid governing 
equation while the liquid governing equation is solved using a 
Control Volume Method (CVM). The liquid-solid behavior is coupled 
and the solution is obtained in matrix form for the whole thermal 


field (solid and liquid) at once. 


1.3 Program Development 


In this section a brief synopsis of the various research activities 


is given in the time frame of the reporting period. 


Period July 1989 - December 1990 Design of the new experimental 


set-up and requisitioning of its major components. Completion of 
the Macor data acquisition with the infrared thermography and 


preparation of the 188-89 final report. 


Period January 1990 - July 1990 Derivation and implementation of 


the CVM method and coupling with the existing BEM for the solid in 


matricial form. Code development and validation with the Macor and 
Aluminum experimental data. Completion of the construction and 
testing of the experimental set-up for the study of dropwise 
evaporation on Macor induced by radiant heat fluxes above the solid 


surface. 


2. EXPERIMENTS WITH RADIANT HEAT FIELD 


2.1 Radiant Set-Up Design Requirements 


The initial portion of the experimental program was developed using 
solid materials heated by conduction from below. The heat source 
was an electric heater placed ender the solid block which was 
simulating a semi-infinite solid with respect to the droplets 
deposited on its upper surface. A more realistic representation of 
a fire environment requires the semi-infinite solid to be heated 
from above by radiant heaters simulating flames generated by the 
combustion of the pyrolysis products released by the solid 


material. 
The design requirements for the experimental set-up are: 


a) uniform radiant heat input from above the 
solid surface capable of maintaining an 
initial solid surface temperature of about 
2004°c 

b) infrared camera unobstructed view of the 
droplet and of the area of surface influenced 
by the evaporative cooling process (up to 5 
droplet diameters) 

Cc) unobstructed access for the droplet depositor 


used in the previous Macor experiments for the 


time required to deposit a single droplet (to 
limit thermal radiation damages to the 
dispenser) y 

a) uniform or linear temperature distributions in 
the solid at the initial steady state to 
satisfy the initial condition required by the 
BEM 


e) radiant heater output control to maintain 


uniform conditions during the test 


These design requirements were carefully considered and a balanced 
compromise has been achieved among the various items. In the 
following these various elements will be discussed and the final 


set-up configuration will be described in detail. 


2.2 Solid Thermal Distribution Code 

The first issue to be considered is the thermal configuration of 
the solid in its initial steady state. A simple computer code is 
developed to predict the temperature distribution in a block of 
Macor subjected to the following boundary conditions: a) radiant 
heat input, convection and re-radiation at the upper surface, b) 
convective losses at the block exposed sides, c)conduction heat 


losses at the bottom surface. 


Note that Macor exhibits a thermal conductivity of 1.3 W/mK which 


is lower than most insulating materials. Therefore, the idea of 
obtaining an adiabatic boundary at the sides of the block is 
unpractical since one would obtain an extended surface 
configuration instead of an insulating barrier. Tot limit. the 
losses by radiation from the solid sides, a gold paint is applied 


which drastically lowers the Macor emissivity. 


The computations performed for a thick Macor slab show that the 
initial temperature profiles in the solid are not linear nor 
uniform. The solution to this problem is achieved by reducing the 
thickness of the Macor and by adding a chill plate at the bottom 
surface to induce a large heat transfer in the vertical direction. 
The radial heat flux distribution is overshadowed by the vertical 
component of the flux and a linear temperature profile is obtained 


in the central portion of the Macor solid. 


2.3 Solid Geometry and Boundary Conditions 


The results illustrated in the previous section dictate that a 
Macor slab should be used. The optimal thickness of the slab has 
been evaluated in the order of one inch. This thickness is large 
enough to retain the semi-infinite solid configuration with respect 
to the droplet and it provides at the same time the maximum area of 


linear initial temperatures required by the BEM. 


The shape of the Macor block is selected as a 6 inches square tile 


(one inch thick) since this configuration is readily available. 
The area influenced by the evaporative cooling process for a single 
droplet is of the order of a circle of 3 inches in radius. 
Therefore, the 6x6 size of the tile provides sufficient surface for 
the experiments. Note that the tile size is limited by the heater 
size and output since a uniform radiant heat condition is 
postulated over the whole tile surface. A reduction in radiant 
heat flux at the edges of the tiles will increase the steady state 
temperature distribution non-linearity. The convective boundary 
condition at the tile edge is the minimum possible heat transfer 
conditions that could be easily implemented. The possibility of 
evacuated spaces at the tile periphery has also been considered but 


is discarded. 


The chill plate configuration does not require addition comments. 
The Macor tile is mounted on the chilled plate with a conductive 
paste to neutralize the contact resistance. The chill plate is 
cooled with tap water which is delivered at a constant temperature 
of 15 °C. The cooling water temperature rise in the chill plate is 
less than one °C. Therefore, the assumption of constant 


temperature at the Macor tile lower surface is reasonable. 


2.4 Heater Design, Control and Configuration 


Three conical heaters are used to provide the radiant heat input to 


the Macor tile. The heaters are the cone calorimeter developed at 


CFR. One cone calorimeter heater was available and the other two 
heaters are similar in dimension and power output but much simpler 


in construction. 


This three heaters system is connected to a three phase 220 volts 
power supply in a delta configuration. Two heater elements are 
identical and the third is very close to the other two. This 


provides a balanced electrical load. 


The power supply is driven by an automatic controller which senses 
the heaters average temperature through three thermocouples and 
reference it to the set-point. This system is able to maintain 
constant heater temperature which results ina stable steady state 
Macor surface temperature. The Macor surface temperature that can 
be achieved with chilled plate cooling water at about 10 °C is in 
the order of 170 °C which is satisfactory for the evaporative 
cooling experiments. The heater temperature required for this 


Macor surface temperature is in the order of 750 °c. 


Higher surface temperatures are achievable if the chilled plate 
under the Macor tile is not in operation. This is at the expenses 
of the uniformity of the initial surface temperature which is a 
condition that must be met for the exact evaluation of the 
evaporative cooling when the experimental data is compared with the 
prediction of the coupled solid-liquid model. This model is not 


used for the nucleate boiling regime of vaporization which occurs 


at temperatures in excess of 160-170 °C. 


2.5 Infrared Image Shielding From Background 


Information on the transient thermal behavior of the Macor tile, of 
the droplet vaporization and on its size are gathered via an 


infrared camera which sends the information to a tape recorder. 


The infrared thermography of the irradiated surface is complex 
since reflection is also present (the Macor emissivity was found to 
be equal to 0.94). To shield the camera input from most of the 
reflective input the portion of the surface of concern is seen 


through a chilled pipe. 


The chilled pipe is constituted of a coiled copper tubing painted 
with a high emissivity paint (approximating a black body) and 
insulated on the outside. This pipe effectively absorbs all the 
stray and diffused radiation thus sending to the camera only the 
direct radiation incoming from the portion of the Macor surface in 


sight. 


With this set-up a steady state uniform temperature profile is 
observed in all the field of view. By repositioning the three 
heaters carefully most of the reflection is kept out of the camera 


yielding meaningful information. 


10 


The transient thermography of the surface Macor surface is 
recorded. At selected time during the droplet evaporation one 
recorded frame is isolated and digitized. This digitized 


information is then used as input to the post-test code. 


2.6 Post Test Code Modifications 

The code previously developed was able to differentiate between 
meaningful information and noise and could provide some curve data 
fitting correlations. 

A new sede is available which doubles the data used in the surface 
temperature curve fitting routine and provides a better resolution 
for the measurement of the droplet size. This code will be used 


for the data reduction. 


2.7 Notes and Procedures 


We are currently studying the effect of placing a water droplet on 
a surface which is heated from three heaters. The radius of 
influence, time of evaporation, and the temperature profile induced 
by the water droplet on the heated surface at various times into 
the evaporation process are all important parameters which will be 
addressed. The setup and the experimental process will be 


discussed. 
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2.7.1 Setup 


The first step in setting up the experiment is to determine a 
proper heater setup. The heaters should produce as uniform as 
possible radiant heat flux on the macor surface. It has been 
GQiscovered that, although the macor has a very high emissivity 
(0.94), there is a non-negligible reflection off of the surface due 
to a very high incident heat flux. When a heater is placed facing 
the infra-red camera, the reflection is increased. Although the 

amount of reflection can be measured to a fair degree of accuracy, 
it is better to minimize the amount of reflection into the camera 
lens. After all the purpose of the camera is to see the emitted 
heat flux, not the reflected heat flux. For this reason, it is 
advisable not to place a heater opposite of the camera. The heaters 
cannot interfere with the placing of the water droplets. The water 
droplets are placed by a device which rolls in and out of the setup 
via two aluminum tracks. The camera is placed across from where the 
droplet dispenser rolls towards the setup. This way, the heaters 
are free to be placed anywhere except across from the camera. Two 
of the heaters are placed facing perpendicular to the camera's line 
OfMvisrqnts This” way reflection is minimized and since they are 
directly opposite of each other, they produce a uniform heat flux 
on the surface. The third heater is placed on the same side of the 
macor block as the camera. It is placed further away then the other 
two heaters. The camera looks through a hollow tube made up of a 
refrigeration coil. The tube is used to allow the camera to look 


directly onto the macor surface while preventing the camera from 
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viewing extraneous radiation from the surroundings. The tube has 
its inside painted black and has cold tap water flowing through it 
in order to minimize the amount of infra-red emitted from its 
surface. The macor block is attached to a plate with cold water 
flowing through it with heat sink compound. This maintains a cold 
temperature underneath the macor block. Ideally, all heat which is 
absorbed by the macor block would be conducted downward into the 
cold plate. If no heat is being conducted across the surface of the 
plate, there is a uniform temperature across its surface. In order 
to aid in this, the sides of the macor block are gold plated, so 
that there is very little heat transferred via radiation from the 


sides. 


Ze/.2 “Procedure 

The first step is to de-gas the water. If gasses are present in the 
water, the gasses will accumulate on the heated surface and effect 
the heat transfer mechanisms taking place during the evaporation 
process. The de-gassing is done by repeatedly freezing and melting 
the water while the air in the water container is being vacuum 
pumped. This needs to be done about once a week. The heaters are 
turned on and their temperature is adjusted so that the surface 
reaches the desired steady state temperature. Obtaining the actual 
surface temperature is done via a procedure described in the 
appendix. The infra-red camera is properly adjusted via calibration 
from a blackbody source. The camera is recalibrated about once a 


week. Once the proper surface temperature is reached and the camera 


V3 


is ready, droplets can be placed. Droplets at sizes of 10,30 and 50 
microliters are placed at various temperatures. About 10 droplets 
of each size at each temperature are placed so that an average can 
be taken. A VCR records the temperature profile about each of the 
droplets as a function of time elapsed. Each of the droplets has 
its temperature profile captured onto a computer disk at several 
different times into the evaporation process. This is done using a 
"Computer Eyes" software package. The runs are then curve fitted 
using a C language program which fits an exponential curve to the 
temperature profile. The program first displays the temperature 
profile vs. location on the macor block. The user then locates the 
position of the droplet using a mouse. The profile then inverts one 
half of the image onto the other half of the image. Since the 
temperature profiles are symmetric about the water droplet, this 
produces the effect of having twice as many data points. Next, the 
user selects two points, between which, a curve is to be fit. The 
curve fitting program produces two coefficients, which are a 
function of screen coordinates. Another program then converts the 
screen coordinate coefficients into coefficients in term of inches 
or millimeters. The evaporation times are measured with a stop 
watch. The radius of influence is obtained from the curve fit to 


the temperature profile. 
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3. MODEL OF DROPLET/SOLID INTERACTIONS 


3.1 Review of Boundary Element Method 


The conduction equation describes the thermal transient in a solid 
and in a motionless liquid. The full description of the problem 
requires the identification of all the boundary conditions and the 
coupling of the liquid and solid regions at the wetted surface 
under the droplet (see Fig.1). The conduction equation can be 
written for both the liquid and the solid, keeping in mind that the 


thermal diffusivity should be changed accordingly, that is: 


o =aV?T (cle) 


The liquid region is bounded by a liquid-vapor interface where the 
mass transfer is coupled with the heat transfer as extensively 
described in the previous reports. This boundary condition can be 


summarized as: 


D 


kos Vale Dol T=) = 0.624, 2] (3.1.2) 


a 


The liquid-solid boundary conditions express the continuity of the 
temperature and the conservation of energy across the interface, 
namely: 


k, VT, = k,VT (Bees 


Psy Grige) 


LD 


The boundary condition at the solid-air interface accounts for the 


convective and radiative heat transfer contributions: 


“ko v0 =.h,(T-Tajiroe(l —-1.7) (3.1.5) 


The axis-symmetric nature of the problem grants that, both in 
spherical and cylindrical coordinates, the gradient of T is zero on 
the vertical axis through the origin of the coordinate system. 

The initial condition can be a linear one-dimensional temperature 
distribution in the solid and uniform temperature in the liquid or 


uniform temperatures in both liquid and solid. 


The droplet induced cooling of a semi-infinite body has a dual 
time-spatial parametric description. One could select as 
characteristic length some length related to the droplet wetting 
surface or the ratio of the solid thermal conductivity to the 
convective heat transfer coefficient solid-to-air or even some 
mixed defined length, such as the so called 'penetration length' 
(1 = [a t]°-?) where the time could be the droplet evaporation time 
and a is the solid thermal diffusivity. Similarly the time scale 
could be the evaporation time or a grouping of solid related 
properties such as [(pe, k) /h?] or again a mixed definition such as 
the square of some length related to the droplet wetting surface 
over the solid thermal diffusivity. The, point of all “thas 
discussion been that non-dimensionalization of the governing 


equation is not univocal and it should not be pursued at this 
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preliminary stage. 


The difficulties encountered in the application of differential 
methods to the solution of the coupled solid-liquid thermal 
transient suggest that a different approach should be considered. 
Two major considerations lead to the selection of the appropriate 
methodology: a) the relevant event are taking place on the surface 
of the domain; b) the nature of the problem with its’ sharp 
localized changes in the thermal gradient requires an integral 


approach. 


The idea of reducing the problem to a surface problem is appealing 
because only the points on the solid (and liquid) surface can be 
monitored experimentally. An infrared thermographic technique will 
be used to monitor the surface temperature during the droplet 
evaporation transient. This technique is non-intrusive and is 
highly suitable to this particular problem. 

Therefore, the computation limited to the surface points is more 
efficient and allows a more precise definition of the noding in the 


region of concern (e.g. at the droplet outer edge). 


An integral methodology is also desirable because the solution at 
a given point (e.g. point temperature) is obtained by superimposing 
a large number of contributions from the neighboring points hence 
the localized, drastic thermal changes of this cooling process are 


smoothed out. When a finite difference technique is used the sharp 


ed 


gradients are locally amplified thus causing the observed 


instabilities in the solution. 
To obtain the desired result, it is necessary to introduce an 
adjoint equation in terms of the Green's function 'G', that is: 


— =-aV’G (30ar 6) 


By multiplying Eq. (3.1.1) by G and Eq. (3.1.5) by T respectively and 
integrating over the domain making use of Gauss' theorem, one 


obtains: 


ff2ge siete Mas Michie tes do weeny eaten 


At this point the left-hand-side is expanded into two volume 
integrals at time t and at time zero respectively. By G at time t 
and on the boundary as the Dirac function, the volume integral at 
time t reduces to T. Further, one can introduce a new variable u 


in lieu of the temperature which is zero in all the volume at the 


initial time (t = 0). For instance one can define u as: 
Ups ian ee (3 eed 
k, 


With these two modifications the final result is achieved in the 


form: 
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u = af ¥$ (GVu-uVG) -Adsdt (een. 9} 


ts 
This integral is performed over a closed surface. However, the 
surface bounding the solid region below the plane z=0 is far ahead 
of the advancing thermal wave. Therefore, the contribution of the 
points on the surface bounding the domain in the far field is 
negligible. With the implementation of the cylindrical coordinate 
system, the equation is further simplified since the angular 
integration can be expressed in terms of a Bessel function. The 
final forms results in a double integral in time and along the 


radius 'r' which can be cast as follows: 


-+ Visi at fa eae 
ict le. — fala lerorn bb, 2g tg * Df 2) exe | araat 
J40Na 00 4at, 4a, 
Crile Woy) 


In order to simplify the task of handling these complex surface 
integrals, a decomposition in two portions is proposed: the value 
of the forcing function or of the unknown function (Vu) is assumed 
to be constant in the small intervals dr, and dt,. This introduces 
formal errors of the order (dr,/R) * and (dt,/t)* into the analysis. 
The various surface integrals can be recast in the following useful 


form: 


where 'W' is a weight matrix and 'f' is the vector of the forcing 


and unknown functions. 


Ag) 


n 
u= WF +Wy° fy (3 lee 


i= 
The summation term is known since it involves previously calculated 


parameters. The second term on the right-hand-side of Eq. (3.1.10) 


contains unknown functions and the specified boundary conditions. 


In general, each solution at a given point depends on the effects 
of all the other surface points. So far the effect of points 
belonging to the same surface have been considered. fThe inter- 
surface dependance will require special treatment as described 
later with the presentation of the final solid-liquid coupled 
model. More immediate consideration must be given to peculiar 
situations: the effect of the forcing function at the point of 
concern and the singular behavior of points at the origin of the 
spacial and temporal coordinates. 

The first particular case of interest is when the point of concern 
coincides with the source point. This corresponds to the maximum 
influence of a source point because the distance is minimal. 

In order to treat this singularity in detail, an analytical 
integration of the general expressions of 'G' and 'VG' is performed 
over the elementary portion of surface represented by one node 
(e.g., (x - Ax/2) < x < (x + Ax/2) ). For the disk and the plane 


the result is: 
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W (rmr) ty) ser ty uf 22 | vx Ste) a, (3.1.12) 

Note that the term At, is consistent with the formulation of 
Eq. (3.1.11) where the forcing function f is multiplied by the 
weight W. By comparing Eq.(3.1.11) with the general formulations 
for u, it is clear that a temporal and spacial integration must be 
performed on the weighing function. Such integration in this case 
is performed analytically in the space and numerically in time. 
The units of the weight W are those of length as expected. 

About the points at time t = t,, these singularities are addressed 
by staggering the discretization of the domain. The time 
discretization is not carried to the time origin (eg. the present 
time ort, = 0). Since®the weight function is characterized by the 
presence of a high peak (Similar to a Dirac Delta) near the time 
origin, the integration over the first time step has to be 
performed with the maximum possible accuracy, and this was done by 
using the most advanced mathematical routines by IMSL. The 
discretization of the integrand must have a resolution such that 
the peak of the function can be properly captured. Consequently, 
the behavior of the weight function for various governing 
parameters versus the recollection time becomes essential. at 
order to describe the weight function for combinations of the 
parameters (Ar,, At,, a), a non dimensional quantity is identified 
as P= Bare ACale The physical properties of different materials 


determine the value of the parameter P, and consequently the height 


al 


of the initial peak in the temporal behavior of the weight 
function. This information can be used to set properly the duration 


of a single time step. 


In light of the previous observations it is important at this point 
to gain some insight in the concept of recollection time. The real 
time is elapsing as the evaporative cooling process takes place. 
The recollection time is zero at the present time and stretches its 
positive axis in the past. Further, the effect of past forcing 
functions fades as real time goes by, hence the weights of the past 
forcing function keep decreasing. Depending on the type of 
material and on the extent of the cooling, the recollection 
(memory) of the system is longer or shorter but in any case limited 


only to a finite number of previous time steps. 


By rearranging the equations for all the points, one can recast the 


problem in the form 
aes oe (3',1 Say 


where A is the modified weight matrix at the present time. This 
matrix is composed of all the weights associated with the unknown 
fluxes and with portion of the identity matrix associated with the 
unknown temperatures. The vector B is a vector that encompasses 
known quantities at present and past times as shown in the 
summation term of Eq.(3.1.11). Further, B also includes the known 
weight-flux products or the known temperatures at the exposed 
surfaces (cap and plane) which are described by the boundary 
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conditions. 


Matrix Structure. Intuitively, the weight function must have its 
highest value when the source coincides with the point of concern 
and should decay as the source moves further away from the point. 
This behavior is confirmed by the numerical values of several test 
cases. The weights can be arranged in a matrix form where each 
column is composed of all the weights associated with the sources 
affecting one single point. Similarly a row is the collection of 
all the weights associated to one single source affecting the 
various points of the domain. This weight matrix is diagonally 
dominant and decays on both sides of the main diagonal as shown in 
Table 3. Further, it is important to note that, as time elapses 
the weight matrix values decay rapidly and in fact after a few time 
steps they become negligible. This confirms that the system 
recollection time extends in the past for a limited number of time 
steps. 

After the preliminary characterization, the weight function must be 
tested in the full range in matrix form to establish its accuracy. 
Some relevant niecerences between the shape of the matrix in the 
case of aluminum and the corresponding one for the case of Macor 


will be outlined in the following pages. 
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3.2 Liquid Droplet Nodalization 


As the temperature changes, the volume of the droplet decreases 
with time because of the evaporation. Under the assumption that 
the radius (R) of the wet area, which is the interface between the 
droplet and the solid, remains unchanged during the entire 
transient, the thickness (a) of the droplet must decrease until the 
droplet is completely evaporated. Another assumption to be made is 
that the droplet has a shape which is a spherical segment (see 


Fig.2), so that the volume of the droplet can be expressed by 
v=2 a (3R? +a?) (seo Ey 


and, knowing V and R, a can be found as 


3 
a= 3V, Per ae feta ened Sr 
~ ir 3 (2. ae 


In cylindrical coordinates - r,z - the nodal points in the droplet 


are determined by 


1c 
Z*+(r-g,R)? = | Zi) e( 2-948) (3.20 
9; 
at i Am! 
Tey 2 tay 
eA £,5Y i ea é £,Y i - (3 32:8) 
4 2 


24 


where 


a 
== ovals 
y= (3.2.5) 
£ == | 500 O<f,<1 (ananG) 
I3= = | 200 O<gjs1 (Bsi2%7)) 


While f, and g; go through from 0 to l) separately, eq. 3.2.3 and eq. 
3.2.4 give the coordinates for all points of the droplet when 
i=2,IM and j=2,JM+1 (see Fig.3). 


Byesoiving Cor r°and@z irom eq: 37203 ands+eqri3.274, tonevobtains 


os a = 
2 Sea ee (3.2.8) 


= 2a, 
and 
43; a B, + B.27;,; (3.2.9) 
where 
_ itg, +2(1-9,) 
Bis = maga? £ (Se2510) 
£,Y iY 
and 
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2| 9;+ 
pice ieee (3.29%a) 
1 = . 
f iy Pil 
a, =1+f,? (3.240 
l-g. 
a, = 218382 - (93+ el (<0) (3.2.13) 
Gaga Bt a ed Loo) Re (3.26 


For points of j=1, i=1,IM+1 which are symmetric to z-axis, 


reteties -r{i,3) (3.20358 


Zli,1)4= ZL (3.2.16) 


For fake points underneath the interface, which is symmetric to r- 


axis, 


zy (IM+1, 75) = ciIM-1y/7) (3. 200 


Z(IM+1,j) = -Z(IM-1,j7) (3.2.18) 


Above the droplet surface there is another layer, i=1 and j=1,JM+1. 
The expression of r;;, 2; can be given by letting f.=1 and 


substituting a, to a in eq.3.2.5, where 
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aj.= a+ (Z(2,2)-2(3/2)] (322 49) 


Or e€q. 3.2.3 and 60, j62 a4 


3.3 Control Volume Method Formulation 
The conduction equation for the droplet without the consideration 


of radiation effects is expressed by the Fourier's law 


Pe aver (es ered) 


where T is defined as the temperature difference between the 


droplet and air: 


ig= Taropiet ~ Tair (3.3.2) 


By integrating over the volume of the droplet and applying the 


Gauss' theorem for the right hand side term of the eq. 3.3.1 


OT si 
ae = ® 
ae” Mee re (3.3.3) 


Noting in cylindrical coordinate, 


dV=2nr da; dS=2nr dl (3.3.4) 


one has 
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oT u ae 
Joe 08 ean (3°35.5) 


L 


or 
r,A2 aa=a see oe CA iE gE WELO ys Fstic, soo Icotrp Touts sa Iptr, 
At AX, 2 AXy 2 AX, 2 Ase 2 


(3.3.6) 

Referring to Fig.4, 
ran MESURE Le Ca a, ae 
Ine = (3.g0eg 
lop = ¥ (53,5 SORE ps ‘ (Zs Zier) (3 eae 
La YATRA; ~en te, De & Zits “#4, FAY (3.3.10) 
A= "07 25 at 2) see atl ae) (3.32 
To= Owoo re urhd, + earns or fei) (3 3: az} 
Ty = ON ZO rt La eee ee eges to) C3: AS ay 
Dy O25 (Lae nee tal pen es ela (3.3.14) 
Teq= On 25.( 25.1, seats ey tie + aes oe) (3.53 Le) 
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Te = 0.25 (Ls, 542 + Zi, 502 * Lier jer * Lier, 542) (3.3.16) 


Ze =.0. 25 Zig, hy tre a eae + Zi, 5) (353847) 

2y = 0, 2502 rn Pee Hey, 7 + Zz 5) (3 aake) 

Zy = 0. AGN Ape 2 et 2.925 + Z,.3: 5) (3.3.19) 

Bes > eT Ugo PS ai fa eV a PO (sme o20) 

Sy 0120 Bee + Let Lies gst 2 ie3, j02) (Smo 20) 

and 

Ro inte) eee (3/3 .22) 

ren aA Wah ena Vk A Oa van ie (63% to 231) 


AX, = SCT rape Zar.) (sito ee 4") 


AZ aa tele il aensp) (an. 25) 


V aNd bea 
Po = 
Ae wedrete gneiss lia, 4 * Sz, 5-17 3, 9-1 + Sir. y2ae1,3 * Fi, 97 i,9 
(Fase 36) 


where 


AN, 


Op ap tae) 


Peet ase ae (Buse 


a: oie ape essa (33.208 


ay rjea — 3.3829 
oe I,dA 2Ax, ( ) 
Est my ap OA A ee 3.anoe 
ag a =(ay 4 tanec a aie) (3. 3a 


Note that all aj;;'s are time dependent coefficients since the 


geometry of the droplet, in terms of r,;; and Z;;, vary with time. 


Applying Crank-Nicolson scheme for eq. 3.3.26, and knowing all T;,;'s 


at time step n+l being only unknowns, one can have 


+1 
Lig ie 
At 


fel n+l +1 n+1 +1 n+1 +1 n+1 +1 n+1 +1 
= 0:5 (a7, qs. Dimtout iva Lina gtes,$ lang teageans ley hair 


n n nN n n n n n n nN 
40.5 (as, 541 TG, 561+ @ie1, 9 Tier, 3449, 977, 3+ 43-1, 3 TH-1, 544i, 3-1 Ti, 5-1) 


(3.3.28 
Rearranging the equation above: 
+1 aL +1 +1 fa 
Dis jer Ti, jor * Dior, 3 TH01,9 + Ds, 9TH 7 + Dyer 5Th-2,.9 7 Py jaTiga = Psy 
(3.358a0 
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where 


Dy yey = 0.5 At at 54 (323234) 
Dates OPS A Canty (sac so} 
De, pec NA tidsery (3.3.36) 
fe eae Slo ecg a (oe 7,} 
Dae CD in et Diyas Oye 161, 7) Garo oh) 


n n n n n n Nn n 
Reg = 0.5A0 (aj, 541 Ti, 561+ Gina, 5 Tis1, 3+ @i-1, Ti-1, 5+ 45, 5-1Ti, 5-1) 


41+ O).5 Morar err 3 (S03 95) 


For points near the edge of the droplet, all above expressions will 
be modified due to the shape difference of the elements (see 


Fig.5). The area dA of the triangle ABC is given as 


GAS 0225 \(1,,* leap. (3.3.40) 


Instead of using averaging procedures for the rectangular elements, 
point O, as well as point N and S, is taken as the center of 


gravity of triangle ABC: 


i 


ym oars (255 7- 2ie1,5 + a2 i41,3 7 PT i, 3) (3.3.41) 


Ly ma Zp BOCE 5 Es >4) (3565004 2) 


where 


cyt 


apn ee ae a (3.3.43) 
0.5 (Ls 5*2y, 501) “Lier, 5 


ie) 


b= 2 (333 244) 


; 5 (2542, 5723, 541) aor 
ip tote atl pe ear 


io) 


Referring to eqs. 3 3270— 3.362570) ete pt ae ee eo Ax, can 


be found. 


Again, applying the Crank-Nicolson scheme, the conduction equation 


for the edge points becomes: 


by 5-1Teegna ti ety Teg ei fF Dia Tin oe tiie eee 
where, Similarly to eqs. 3.3.34 - 3.3.39, 
bee, © 0,5 Atal sy (3.3.46) 
by yin -0tS A bay (3.3.47) 
Die eae O oA bas (3.3.48) 
Dynan tl, LOR eee tis ae) (3.3)44g7 


RiG =0,5A8 (iy 9Ti-1,9+@f 7-27 ty t+4in, sTA DO) +(14+0.5Ataz;) T;; 


(3.3.50) 


leotr +r) 
a AB\‘*A “B 
2;.,;,7°= o-oo 3.03 so 
ees Ida 20%; ( ) 
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aj. = > WAS pn eo 


a Sane. ee 7 (37.1315) 


Qpeemes (ayes Fey yet 4i37, 7) C353b 54) 


The governing equations for points on the droplet surface and 
points which are symmetric to Z-axis are determined by boundary 
conditions. 

Due to phase change and convection, there are both heat flux and 
mass flux loss through the liquid surface; therefore, by taking 


energy and mass balances on the surface, one has: 


io Lien | siete My pa 323455 
a c=. |Pabad se + B(Ty-T,) ) 


After simplifying and discretizing, the equation above becomes: 


2 
Te Ty o.c24h A( D3 Ae). p( BaeTis a, (3.3.56) 
2 


Ax K, C,\ a 


a 1 


or in more general form 


Dy Taig thy 344 = Rj CSkiay 
where 
Axh Axh 
Dra ae : Do = l= (on oO) 
set 2K “el 2K, 
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R <= 0.624Axh A/{ D Xi -Xs _ Axh,, (3.3.59) 
1.9 K CHE 1-x Ke 


W a 


Referring to Fig.6, 


Ax = (r--r*)?2+ (Z°-Z*) 2 (3.3.60) 


I=) 0% 25,07, Bt ht Foggy Tonge (343,. 62) 
TH =9032 52, othe gm +p? eyst) (3.3.62) 
ZirstOe2o(Zilwme, See sam + Ze } (328ueuy 
Z yre0 2502, ,t2ob pada t Zia) (3.3.64) 


In addition, because of the symmetry of the temperature field about 


the center line of the droplet, the boundary condition is given by 


OT 
—_ = 3.35 Oo 
or ss 0 ( ) 


which can be further expressed as 
Ds 4T3,1 + 53,2Ti,2 = Ris (3.3.66) 
where 
b;,=1; b 


52-1;  -R; ,=0 (3.3 6a) 


Because of the continuity of the heat flux at Z=0, the boundary 


condition for liquid and solid interface is: 


Sean : ba ee (3.3.68) 


Further discussion will be given in next section. 
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3.4 Coupled Model Description 
In order to couple the solid and liquid parts together, the 


interface between solid and liquid must be focused on. 
From the analysis of the boundary element method on Section 3.1, 


the equation for all solid surface points can be written as 
ug =) fF Wey + Ty (3.4.1) 
where 


sen Dh ID a (3.4.23 


By definition 


q h 
us Se cae = Loy C18 sat bra (ota 03) 
and 
du 
ie =-—_ Belo louht 
az ( ) 


fe a7 h, h, 
is Se (PS 3 CS ay Mal 
AZ x, | s a rahe ‘3 ( ) 


For interfacial points: 


aber Oe, (2 ane) 


Applying boundary condition 3.3.68, the above equation becomes 


mE) 


K, oul s 
= —- f er = 3 © 4 e v 
a ses K, (Ts Ta) 


where K//K, is the ratio between water and solid conductivities. 
Again, letting T to be the temperature difference between solid 


surface and air, and referring to Fig.7, eq.3.4.1 can be written as 


JM K aT h JM+JE h JM+JE 
T, = ralee hy ae T,W,, +2, + (T,-T,)|1+— Wer 
: d, K, aZ k , Ks 2, z : i F K, D : 
(3.4.8) 
Knowing that, for interfacial points: 
© Tyy,5% Tinea, 5 
T, = 5 (3,.459) 
and that, for exposed surface points: 
Tatler (3.4.10) 
and: 
SS = Dyy-1,x~ 11, k (3. 45a 
daZ} , AZ, 
one has the equation for points at the interface: 
Wa: u/s eae W,, K 
0.5-—2 i) Toe PVPS ee (0.54 dd Zz)? 
| AZ, K ay » Kas AZ ee 
JM/k K Ws, JM+JE h JM+JE 
A ee ee K bers Pe at W.T alr +4(T -T )i1+— W. 
IM,k k+IM,k 
D K, AZ, Ms raked ‘eae : ; ‘ K, yD cis 
(32 ie 
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and for points on the exposed surface: 


JM JM 

K, Wy K,, W, h 
VET x dG eP Lig Ai a 
& K Tyy-1,k > K Ag IM,k ( K, 5/3) IM,J 


= s 


en, JM+JE/j JM+JE 


h 
» Welty ke = ryt (tet) 1+ SS A (Same) 


Kg k=JM+1 Ss k=2 


The right hand side terms of eq. 3.4.12 and eq. 3.4.13 are known 
values which contain the results from previous calculations. 

After assigning the governing equations for each point of the grid 
mesh, either from the conduction equation or boundary conditions, 


one finally has: 


Bio Tec= ok (Seas) 


where T is the unknown temperature difference vector at time step 
n+1 which has to be found; B is a matrix containing all the 
assigned coefficients, while R is a right hand side vector which 


includes the known values (see Fig.8). 


By solving the inverse of the matrix B, the solution for T at n+l 


time step is given as 


T = B1l-R (374 5) 


Based on the newly calculated temperature field at current time, a 


new droplet volume can be computed from 


oe 


Ve ee Void mary (3.4.16) 
where 
my 
av _ 2n(0.624)h; D)3 [aera (3.4 Ale 
dt pC v4 aoe 
wa a 0 a 


in which dl is an integration along the droplet surface. 


A flow chart of the procedure for solving this coupled problem is 


given in Fig.9. 
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3.5 Special Issues and Remarks 


3.5.1. Recollection time. 

The recollection time is a parameter that has to be fixed at the 
beginning of each run in relation with the solid surface material. 
The "memory" of aluminum is totally different from that of Macor: 
more specifically, the sensibility of aluminum to the heat fluxes 
from its past history is in the order of magnitude of 1 second, 
while the corresponding "memory" of Macor stretches back for more 
than 10 seconds. Therefore, a recollection time of 10-20 seconds 
can be imposed to comprehend conservatively the behaviors of both 


aluminum and Macor. 


3.5.2. Normalization of the matrix. 

The matrix of the coefficients of the coupled model has to be 
properly normalized in order to avoid any kind of numerical problem 
during the inversion procedure. This is due to the composition 
itself of the matrix: since it is constituted by terms of very 
different kind and order of magnitude, its inversion without a 
preliminary normalization would generate an increase in the 
sensibility to the initial and boundary conditions and a decrease 
in the level of accuracy of the solution. Therefore, since the 
highest values are located on the main diagonal, the matrix was 
normalized by dividing each row by the term belonging to the main 


diagonal, and the same operation was performed on the vector of the 


wo 


known terms. 


3.5.3. Setup of initial conditions for the temperature field. 

Before the beginning of the transient, the liquid temperature is 
imposed to be equal to the ambient temperature (20°C) everywhere, 
including also the fake layer underneath the liquid-solid 
interface. The huge heat flux corresponding to the sudden start of 
the transient creates a discontinuity that is strongly sensed in 
the numerical simulation: in the heat flux behavior, this results 


in an initial damping, which becomes smoother with time. 


3.5.4. Edge of the droplet. 

The largest variations in the temperatures and heat fluxes of the 
liguid occur near the edge of the droplet. In this region, some 
improvements in the numerical simulation could be obtained with a 
more refined integration path, since the nodalization chosen does 
not reproduce the shape of the edge as well as it reproduces the 
central part of the droplet. However, this kind of modification was 
aes eae tried and the corresponding improvements were 
observed to be negligible if compared with the complexity of 


calculations required. 


3.5.5. Heat transfer coefficient. 
The overall heat transfer coefficient used for the liquid surface 
is based on experimental measurements, from which a constant value 


along the whole droplet surface is given. At the same way, also the 
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heat transfer coefficient for the solid surface is assumed to be 
constant. These assumptions are perfectly reasonable for the 
regions where separation between only two surfaces occurs (solid- 
liquid, solid-air and liquid-air), but they create a mathematical 
discontinuity at the intersection points between all the three 
surfaces (at the droplet edge defined above). The sudden transition 
from a region of high fluxes (solid-liquid) to a region of almost 
negligible fluxes (solid-air) seems to behave like a step function 
in the numerical simulation, and this does not correspond to the 
physical reality of the phenomenon. Some fluctuations in the 
temperatures and in the heat fluxes near the edge of the droplet 


can be easily explained from this point of view. 


3.5.6. Sensibility of the solution scheme and optimization of the 
matrix of coefficients. 

The matrix of coefficients in the coupled model is characterized by 
a very particular shape, as it was previously described. This shape 
is univocal determined from a qualitative point of view, but a lot 
of input parameters, that can be chosen with a certain degree of 
freedom, can create some numerical discrepancies between two 
solutions of the same case. This was observed and analyzed in order 
tow -obtainslas-sort., of,joptimization«~of. the,matrix..in. terms of 
combination between accuracy and CPU time. The results of this 
analysis can be summarized as follows: 

1 - The meshes in the liguid grid play different roles in the 


solution scheme: since the heat fluxes are predominantly axial, the 
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vertical nodalization is more important than the horizontal one to 
achieve a better understanding of the temperature distribution 
inside the droplet. 

2 - The points on the solid exposed surface are characterized by 
temperature gradients which decrease with the distance from the 
edge of the droplet; therefore, it is possible and useful to impose 
a discretization that starts with the same level of accuracy chosen 
for the region under the droplet and decreases in proportion with 
the distance r. 

3 - The matrix lines related to the solid part of the model 
contain a larger number of terms than the lines related to the 
liquid part, and the basic characteristics of the involved 
equations show that, in terms of CPU time, an increase of the 
number of solid points is more expensive than the same increase in 
the liguid region. 

4 - The liquid-solid interfacial points are the most delicate 
ones: their behavior, in terms of accuracy of prediction of 
temperatures and fluxes, determines the convergence of the whole 
solution. Their number, however, has to be the same that is chosen 
for the horizontal discretization of liquid and solid separately; 
therefore, it is important to set that number of horizontal nodes 
in such a way that the significant variations of the interfacial 
temperature are individuated without a too fine discretization. 
The nodalization which was obtained as an optimized result is 
characterized by the following Sobdivtwilonett 


- liquid height: 15 layers; 
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- solid-liquid interfacial segment: 6 nodes; 
- exposed solid surface: 24 nodes, with the distance between them 


doubled every 6 nodes. 


3.5.7. Numerical problems from VAX to IBM. 

The code was initially implemented on the VAX3100 of the University 
of Maryland at College Park, and only in a second time it was 
implemented also on the IBM 370 at the same location. This was 
supposed to be a trivial transfer operation between two well known 
machines: in the reality, that operation turned out to be a complex 
procedure involving a certain number of numerical problems not 
always easy to be solved. 

Some mathematical subroutines from the IMSL library, for instance, 
were found to behave differently on the two machines in relation 
with the number of elements of the matrix of the coefficients and 
of the vector of the known terms. The inversion of the matrix was 
the most delicate step, also because the different intrinsic level 
of numerical accuracy from VAX to IBM is per se a cause of little 


discrepancies in the corresponding truncation errors. 


3.5.8. Inconsistency of the theoretical contact temperature. 

The experimental and numerical results obtained for very different 
solid surface materials demonstrate that the theoretical contact 
temperature proposed by Seki et al. is a very rough approximation 


of the real distribution of temperature that characterizes the 
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interfacial behavior during the evaporative transient. The concept 
of contact temperature appears to be useful in terms of stability 
of the solution, but its numerical value is close to the real one 
in a very limited number of cases. 

The demonstrated inadequacy of the theoretical contact temperature 
represents a turning point both in the calculation of the 
interfacial temperature distribution and in the definition of the 
radius of influence for the droplet. An important part of the 
analytical program of the next grant period will consist of 
substituting that contact temperature with a new reference 


parameter based on the interfacial heat flux. 


3.6 Code Validation 


The coupled model presented in 3.4 was used to reproduce 
numerically a large number of experimental tests provided by the 
infrared apparatus described in section 2. A direct comparison 
between numerical predictions and experimental results is possible 
for the evaporation time and for the distributions of temperature 
under the droplet and on the exposed solid surface. 

The evaporation time of droplets of different size (10 ml, 30 ml 
and 50 ml) is the variable which gives the most consistent insight 
into the prediction capability of the code. Therefore, a very 
large number of data from tests on aluminum and Macor was 
collected, in a broad range of initial solid surface temperatures, 
and all the corresponding evaporation times were calculated by the 
code and eventually compared with the experimental values. The 
experimental values reported here are, for each case, the result of 
an average of the slightly different evaporation times obtained 
during more repetitions of the same test: this holds the margin of 
the experimental errors under an acceptable 5 percent. 

The results of this comparison are presented in the following page 
(Table 1), and the remarkable agreement between experimental and 
numerical results is showed in Fig.10 and Fig.11, where the 
calculated evaporation times on aluminum and Macor are plotted 


versus the corresponding experimental values. 
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Table 1 - Evaporation time: comparison between experimental and 


numerical results for aluminum and Macor 
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Aluminum 


Specific Heat (J/Kg°C) 962 
Thermal Conductivity (W/m°C) 180 
Thermal diffusivity (m*/s) 4.55x10” 


Table 2 - Solid Properties 


Macor 


835 


6.19x10°’ 
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The excellent agreement between experimental data and numerical 
predictions is evident, and a point that has to be outlined is the 
coherent degree of accuracy in the numerical simulation of tests 
performed on two materials extremely different from the point of 
view of their physical properties (see Table 2). 

The experimental measurements via infrared camera also show that 
there is a sharp increase in the evaporation rate during the final 
part of the transient, and this is correctly reproduced by the 


numerical simulation. 


The predictions of the distributions of temperature in the droplet 
and on the solid surface constitute a main test in the validation 
of the code, even if the experimental data are obtainable only for 
the solid. The calculated distributions of temperature in the 
liquid, however, are very useful to demonstrate that the axial heat 
flux constitutes more than 90 per cent of the total heat flux 
during the transient. In the tests on Macor, it is possible to 
individuate some points at the edge of the droplet where the axial 
heat flux drops below that 90 per cent. This is to be expected if 
one considers the SO ara of the interfacial temperature 


distribution (see Fig.12). 


About the temperature distribution in the solid, there are two 
basic fields of investigation: the first one consists in the 
analysis of the behavior of temperature in the solid area below the 


droplet, the second one consists in the analysis of the cooling 
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effect of the droplet on the solid "far field". Both the subjects 
have been investigated in this research, and the main results are 


presented here. 


Before the coupled model was obtained, the solid temperature below 
the droplet was a boundary condition and it was imposed to be 
uniformly equal to the theoretical contact temperature proposed by 
Seki. With the coupled model, the matricial equation that 


constitutes the solution procedure was set in the following form: 
AX=B ck eed Sea ot 


where A is the matrix of liquid and solid coefficients, X is the 
vector containing the unknown terms and B is the vector containing 
the known terms. The coupling of the liquid and solid models 
requires that the unknown terms in X are temperatures only (in lieu 
of fluxes under the droplet and temperatures on the exposed 
surface, which were the unknowns in the solid model). With this 
solution scheme, the calculated temperatures under the droplet were 
compared, case by case, with the corresponding theoretical contact 
temperatures: the results showed the inadequacy of the concept of 
contact temperature, both when the solid surface consists of a 
material with very low thermal conductivity like Macor (that is the 
lowest conductivity case considered here: k=1.29 W/m C) and when it 
consists of an “opposite" material like aluminum (k=180.0 W/m C). 
It is particularly significant to compare the temperature 
distributions obtained for "extreme" materials like aluminum and 
Macor in cases when the theoretical contact temperature is the 
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same. One of these comparisons is showed in Fig.13. 

The temperature distribution on the exposed solid surface is the 
parameter that allows to estimate how far the cooling effect of the 
droplet is sensed by the solid. To locate the boundaries of this 
area of influence, it is necessary to define a reference 
temperature, T..,, that has to be compared with the current solid 
temperature T during the transient: of course, it is totally 
arbitrary to set the limit, in the difference (T-T,,,), under which 
the corresponding points of the exposed surface can be considered 
not influenced by the cooling effect of the droplet. 

At the beginning of these researches, it was decided to define a 
"radius of influence" by using the theoretical contact temperature 
as a reference onpererane: More precisely, given the initial 
solid temperature T, and the theoretical contact temperature The 


was defined "radius of influence" r, the distance from the axis of 


symmetry where (see also Fig.14): 


Ta CeO el Ten (2. 6ea. 


This definition was used to plot the calculated radius of influence 
versus the experimental value, and the results for a typical test 
on Macor are presented in Fig.15. 

The above mentioned inadequacy of the concept of contact 
temperature affects also the definition of the radius of influence, 
but this is more evident when the solid surface consists of a 
material with very high thermal conductivity. In fact, the 
distribution of temperature in the experiments on Macor is always 
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characterized by sharp differences between the values under the 
droplet and the values on the exposed surface, and also by strong 
gradients aT/dr on the exposed surface: with this configuration, 
any reference temperature is apparently acceptable in order to 
individuate an area influenced by the cooling effect of the 
impinged droplet. This is not acceptable, of course, from the 
point of view of the physical meaning. Moreover, the same method 
becomes highly inaccurate when one tries to apply it to solid 
surfaces like aluminum, since the differences of temperature and 
the radial gradients, in those cases, are far below the minimum 
experimental error (At=0.1-0.5 °C versus exp.error=1.0-2.0 °C). 
To find a correlation between a radius of influence and a 
meaningful physical parameter, it was eventually decided to use the 
fluxes in lieu of the temperatures. Since it is: 

(ease) 2 pelcky (ae 

A,t R or 

where A, is the wetted area, it is also possible to find the points 
where, for example: 


i pVh,zR edn 


0.0 “a 
A,tk or 


(3.6.4) 


On this conceptual basis, a new definition of the radius of 
influence is currently being developed. The use of heat fluxes 
instead of temperatures as reference parameters is expected to be 
much more effective to obtain a more general description of the 


cooling effect of the droplet on various materials. 


a. 


A simplified model for droplets evaporating on high conductivity 
solids was based on the assumption of uniform and constant 
temperature under the droplet and on the dominance of the axial 
component of the flux over its radial component in the liquid 
droplet. The results for Macor, however, show that the assumptions 
of the simplified model cannot be applied to low conductivity 
solids since the solid-liquid interfacial temperature is neither 
constant nor uniform throughout the evaporation process. In Fig.16 
the solid-liquid interfacial temperature is shown for different 
times and locations under the droplet. The results presented refer 
to a 30 pul droplet which, deposited both on aluminum and Macor, 
generated the same calculated contact temperature of 82 °C. It is 
clear that the hypothesis of uniform and constant temperature at 
the solid-liquid interface is acceptable for aluminum but not for 
Macor. 

As it was previously anticipated, the assumption of one dimensional 
heat flux in the droplet axial direction (normal to the solid 
surface) can be considered quite accurate (already shown in 
Fig2i2,)% 

A set of plots obtained with the current version of the code is 
presented in Figs.17-20. In these figures the surface temperature 
distributions and the interfacial heat fluxes for aluminum and 
Macor are shown. They constitute the most recent and reliable 
numerical results obtainable with the coupled model described in 


this report. 


52 


4. CONCLUSION AND FUTURE WORK 


4.1 Remarks 

The coupled model which was obtained in the last year of researches 
can be seen both as a basic result and as a turning point in the 
. numerical part of this work. Even if some characteristics of the 
code have still to be improved or modified, the numerical 
predictions concerning the evaporation time and the temperature 
distributions are reliable and can be used for a broad range of 
materials constituting the solid hot surface. 

With the new definition of the radius of influence, also the 
cooling effect of the droplet will be known in detail regardless of 
the chosen solid surface. 

Therefore, all the results obtained in the studies of a single 
droplet behavior are ready to be averaged and used to analyze the 
evaporation of a more complicated system, consisting of a larger 
number of droplets randomly impinging on a defined portion of a hot 
solid surface. The experimental and theoretical programs related 
to this subjects will constitute the main work of next grant 
period, together with the work that is still required to improve 
the accuracy and the generality of the already existing numerical 
code. These programs of future work are described in the following 


paragraphs. 
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4.2 Experimental and Analytical Programs 


The coupled solution for the low thermal conductivity solid case 
will be improved and compared with experimental measurements of the 
solid surface temperature via infrared thermography. A new set of 
experiments on a quartz semi-infinite solid are under way and 
initial efforts to set-up the infrared thermography apparatus have 
been made. Another area of research that will be developed in the 
next grant period is involved with the generation of small droplets 
(order of one microliter). These droplet are in the same size 


range of those generated by commercial sprinklers. 


The solids analyzed in this initial portion of the research program 
are heated from below by an electrical hot plate. In the actual 
fire application the heat is transmitted to the water droplet and 
to the surface by radiation from above. The new experimental set- 
up where both the droplet and the surface will receive radiative 
heat from gas fired panels will be the main experimental apparatus 
in the next grant period. A solid block with embedded 
thermocouples will also be tested in this radiant heat set-up to 
validate the modelling of the solid thermal behavior for internal 
points. 

The experimental program on radiative heating will be conducted in 
parallel with the development of a modified numerical code; this 
will take into consideration the radiative contribution as a new 


boundary condition replacing the conductive one. 
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At the same time, the results of the experimental and numerical 
work done will be newly used in order to provide a more extensive 
model which will analyze a multi-droplet evaporative transient. 
The data obtained for the single droplet will be averaged and 
applied to the more complex situation that has to be considered 
when a hot solid surface is wetted under the action of a multi- 
droplet dispenser. The project and the experimental tests 


concerning that dispenser are already in progress. 


aie. 
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Fig.10 Evaporation time on aluminum: computations versus 
measurements 
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Ficg..aL2 Axial and radial heat flux components: validation of the 
simplified thermal conductivity model assumptions 
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MACOR (Ti = 119 °C) 


SOLID SURFACE TEMPERATURE (°C) 


isda We Pa Solid surface temperature profiles for droplets on Macor 
and aluminum with identical calculated contact 


temperature 
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Fig.16 Solid-liquid interfacial temperature: simplified high 
thermal conductivity model assumptions 
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C.....-PROGRAM FOR SOLVING A CONDUCTION HEAT TRANSFER PROBLEM IN A LIQUID 
C.....e-DROPLET LAID ON THE SURFACE OF AN INFINITE LARGE SOLID. 
C......-THE INTEGRATED CONTROL VOLUME METHOD IS USED IN THIS TRANSIENT 
C......PROBLEM. 
C......CRANK-NICHOLSON SCHEME IS USED 
PARAMETER (IM=16,JMe1l1,JS=JM-1,JE=4*JS,KNT=10,NIMSL=IM*JM+JE) 
DIMENSION TN1(IM*JM+JE) ,TR(IM*JM+JE) 
DIMENSION TX(IM*JM+JE),TA(IM*JM,4),TB(IM* JM+JE, IM*JM+JE ) 
DIMENSION W(JM+JE-1,JM+JE-1),FTIM( JIM+JE-1,KNT) ,F(JM+JE-1) 
DIMENSION RHS(JM+JE-1) ,CT(JM) ,WO(JM+JE-1,JM+JE-1) 
DIMENSION EWIS(JM),SX(JM),SX2(JM) 
DIMENSION R(IM+1,JM+1),Z2(IM+1,JM+1) ,DZ(JM) ,RSO( JM+JE) 
DIMENSION RO3(IM),2Z03(IM) ,NPTT(50) 
DIMENSION IWK(NIMSL) ,WK(NIMSL+NIMSL*(NIMSL-1)/2) 
COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NDP,NP,TN 
COMMON /INTG/ ATO, ERRREL,IRULE 
COMMON /SOLID/ SK,SC,SR 
COMMON /WATER/ COEFKW,CPW,RHOW 
PI=4.*ATAN(1.0) 
IRULE=2 
AT0=#0.001 
ERRREL=0.001 
NTCI=l1 
NTCS=1 
DTIME=1. 
NTIME=500 
READ(5,*)VO 
READ(5,*)BETA 
RR=BETA*EXP(1./3.*LOG( .75*VO/PI) ) 
UU@EXP(1./3.*LOG(3.*VO/PI+SQRT(RR**64+9./PI**2*VO**2) ) ) 
AA=UU-RR**2/UU 
C..-AIR TEMP. 
TMPAIR=20. 
C..-INITIAL TEMPERATURE OF THE LIQUID DROP 
TEMP0=20. 
C...-INITIAL TEMPERATURE OF THE SOLID 
READ(5,*)TSLDO 
C...-ALL TEMPERATURES USED HERE, TN1l AND TX, ARE DELTAT : U=TEMP-TMPAIR 
TAIR=TMPAIR-TMPAIR 
C..-SOLID CONDUCTIVITY (SK:W/MM K; SR:KG/MM°3; SC:J/KG K) 
SK=180.E-3 
SR=2771.E-9 
SC=962.32 
RHK=SH/SK 
C...DIFFUSIVITY OF SOLID (MM°2/S) 
ALFA=SK/(SR*SC) 
C..-WATER DENSITY (KG/MM“3) 
RHOW=997 .6E-9 
C...WATER CONDUCTIVITY (W/MM K) 
COEFKW=.613E-3 
C...-WATER SPECIFIC HEAT(J/KG K) 
CPW=4179. 
C..-DIFFUSIVITY OF WATER: ALPHA=KW/(RHO*CP) (MM“2/S) 
ALPHA=COEFKW/( RHOW*CPW ) 
C...-AIR SPECIFIC HEAT (J/KG K) 
CA=1009. 
C..-AIR PRESSURE (BAR) 
PA=1.014 
C....-CONTACT TEMPERATURE 


fhe, 


CALL CTEMP0(1IM,JM,JE,TMPAIR,TN1,TSLDO,CT) 


Cee ee cas enc ctetetclenane cietele sc cre crele ausrere, ofevatestotale stederalarerererencheratetenereners snsiatene 


CALL RSD(JM,JS,RR,RSO) 
TTT=.5 

NDP=JM-1 

NP=JM+JE-1 

TNeKNT 

DTO=DTIME 
DRO@=RR/(1.*JM-1. ) 


Co clcrera acter svetetcheterstoncre elele. cre sl el state erate epenslele ele clefereterslelalertersterereretetclerenstete 


CALL WEIGHT(TTT,WO,RSO) 


(oy AS A Aloe a S.8 GOA Gia ame 


DO 1 NT#1,NTIME 


Cut « cles olcne exele GheveLoLals (ce cleloteteteneto eletetelelolelerslcletateteretererelatenerersra lenens! seta tsretete 


TIME=NT*DTIME-.5*DTIME 


C...SET NEW GRID POINTS BASED ON NEW RR,AA 


Cres 


90 


Crorereie 


Cloterote 


Carers 


CALL RZ90(IM,JM,RR,AA,R,Z) 
CALL TRIO(IM,JM,R,2,RO3,203) 
CALL DZ0(IM,JM,R,2Z,RO3,203,D2) 


DO 90 I=1,IM*JM+JE 
DO 90 J=a1l,IM* JM+JE 
TB(I,J)=0. 

CONTINUE 
DO 100 Ie#2,IM-1 
DO 100 J=2,JM-1 
CALL NEWSO(4,IM,JM,R,Z,RO3,2Z03, 

I,J,RN,2ZN,RE,ZE,RW,Z2W,RS,2ZS,RO, ZO) 

ALAB=SQRT((R(I+1,J4+1)—-R(I,5+1) )**24+(2Z(I4+1,3+1)-2(1,54+1) )**2) 
ALBC#=SQRT((R(I,J+1)—-R(I,J) )**24+(Z2(1I1,534+1)-Z(1I,9) )**2) 
ALCD#SQRT((R(I,J)—-R(I+1,J) )**24+(Z(I,5)-Z(I4+1,J) )**2) 
ALDA=SQRT((R(I+1,J9)—-R(I+1,J41) )**24+(Z2(I14+1,59)—-Z2(14+1,5+1) )**2) 
AREA=.25*(ALAB+ALCD) * (ALBC+ALDA) 
DXN=SQRT( (RN-RO) **2+( ZN—ZO) **2 ) 
DXW=-SQRT( (RO-RW) **2+( ZO—ZW) **2) 
DXS=-SQRT( (RO-RS) **2+(ZO-ZS) **2) 
DXE=SQRT( (RE-—RO) **2+( ZE-—ZO) **2) 


‘ IJ=(I-1) *IM+J 


CONST=.25*DTIME*ALPHA/(RO*AREA) 
AIM1J=CONST*ALBC*(R(I,J+1)+R(I,J))/DXN 
AIJP1=CONST*ALAB*(R(I+1,J3+1)+R(I,J+1) )/DXE 
AIJM1=-CONST*ALCD*(R(I,J)+R(I+1,J))/DXW 
AIP1J=-CONST*ALDA*(R(I+1,J)+R(I+1,J+1))/DXS 
AIJ=-(AIM1J+AIJM1+AIJP1+AIP1J) 


TB(IJ,IJ+1)=-AIJP1 
TB(IJ,IJ+JM)=-AIP1J 
TB(IJ,IJ-JM)=-AIM1J 
TB(IJ,IJ-1)=-AIJM1 
TB(IJ,IJ)=1.-AlJ 


IF(NT.EQ.1) THEN 
TA(IJ,1)=AIJP1 

TA(IJ,2)=AIP1J 

TA(IJ,3)=AIM1J 

TA(IJ,4)=AIJM1 

ENDIF 
TR(IJ)@=TA(IJ,1)*TN1(IJ+1)+TA(IJ,2)*TN1(IJ+JM) 


$+TA(IJ,3)*TN1(IJ-JM)+TA(IJ,4)*TN1(IJ-1) 
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$+(1.-(TA(IJ,3)+TA(IJ,4)+TA(IJ,1)+TA(II,2)))*TN1(IJ) 
C..SAVE FOR NEXT TIME STEP 
TA(IJ,1)=“AIJP1 
TA(IJ,2)=AIP1J 
TA(IJ,3)#AIM1J 
TA(IJ,4)#AIJM1 
100 CONTINUE 
C...-FOR J=JM 
J=oJM 
pO 110 I=#2,IM-1 
CALL NEWSO(3,IM,JM,R,2,RO3,2Z03, 
$ I,J,RN,ZN,RE,ZE,RW,2ZW,RS,ZS,RO,2Z0) 
ALCA=SQRT((R(I,J+1)—-R(I+1,J5) )**24+(2(1,3+1)-2(I+1,5))**2) 
ALAB#=SQRT((R(I,J+1)—-R(1I,J) )**24+(Z(1I,J+1)-2(1,J) )**2) 
ALBC#SQRT((R(I+1,J)—-R(I,J) )**2+(Z(I41,J3)-Z(1,5) )**2) 
AREA=.25*(ALCA+ALAB) * (ALBC+ALAA) 
DXN=SQRT( ( RN-RO) **2+( ZN—ZO) **2) 
DXW=-SQRT( (RO-RW) **2+(ZO-ZW) **2) 
DXS@-SQRT( (RO-RS) **2+(Z0-ZS) **2) 
IJ=(I-1) *IM+JI 
Cisse 
CONST=.25*DTIME*ALPHA/(RO*AREA) 
AIJP1=0. 
AIP1LJ=-CONST*ALCA*(R(I+1,5)+R(1I,J3+1) )/DXS 
AIM1J=CONST*ALAB*(R(I,J+1)4+R(1I,J))/DXN 
AIJM1=-CONST*ALBC*(R(1I,J)+R(I+1,J9))/DXW 
AlJ=—(AIM1LJ+AIJM1+AIJP1+AIP1J) 
Cores 
TB(IJ,IJ+1)=-AIJP1 
TB(IJ,IJ+JM)=-AIP1J 
TB(IJ,IJ-JM)=-AIM1J 
TB(IJ,IJ-1)=-AIJM1 
TB(IJ,IJ)#1.-AlJ 
C TB(IJ,IJ)@1.-(TB(IJ,IJ-IM)+TB(IJT,IJ—-1)+TB(IJ,I1IJ+IM)+TB(IJ,I1I+1) ) 
IF(NT.EQ.1) THEN 
TA(IJ,1)=AIJP1 
TA(IJ,2)=#AIP1J 
TA(IJ,3)=#AIM1J 
TA(IJ,4)@AIJM1 
ENDIF 
TR(IJ)@TA(IJ,2)*TN1(IJ+JM) 
$ _ #TA(IJ,3)*TN1(IJ-JM)+TA(IJ,4)*TN1(IJ-1) 
$ +(1.-(TA(IJ,3)+TA(IJ,4)+TA(IJ,1)+TA(IJ,2)))*TN1(IJ) 
C..SAVE FOR NEXT TIME STEP 
TA(IJ,1)*#AIJP1 
TA(IJ,2)#AIP1lJ 
TA(IJ,3)2#AIM1J 
TA(IJ,4)eAIJM1 


110 CONTINUE 


I=l 
DO 120 J=2,JM 
TJ=.5*(TN1(J)+TN1(J+JM) )+TMPAIR 
OSPR beg 
IJ=J 
IF(J.EQ.JM) THEN 
RM=RO3(I) ~ 
ZM=203(1) 
RP=RO3(I+1) 
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120 
C..FOR 


ae RS 


136 


a Ae | 


130 


ZP=2Z03(I+1) 
ELSE 
RM=.25*(R(I,J)+R(I,3+1)4+R(I+1,34+1)4+R(I+1,J9) ) 
2M=.25%*(Z2(1,09)4Z2(1,0+1)4+2(I4+1,3+1)4+Z(1+1,J) ) 
RP=#=.25*(R(I+1,J)+R(I+1,J34+1)+R(I+2,3+1)+R(I+2,7) ) 
ZP=.25*(Z(141,59)4Z2(I141,J341)42(14+2,541)+2(I4+2,3)) 
ENDIF 
DZ1=SQRT( (RM-RP) **2+(ZM—ZP) **2) 
TB(IJ,IJ+JM)=1.—.5*COEFFH*DZ1/COEFKW 
TB(IJ,IJ)#-1.-—.5*COEFFH*DZ1/COEFKW 
CONST1=.624*DZ1*COEFFH*ALAMDA*EWIS(J)*SX(J)/(COEFKW*CA) 
TR(IJ)=CONST1—-TAIR*DZ1*COEFFH/COEFKW 
CONTINUE 
I=IM 
I=IM 
DO 130 J=2,JM 
IJ=(I-1)*JIM+J7 
IF(NT.GT.NTCI ) THEN 
DO 135 JJ=2,J5M 
IJl@(I-1)*JM+JJI 
IJ2=(I-2)* JM+JJ 
IF(IJ1.NE.IJ) THEN 
TB(IJ,1J1)=RKLS*W0(J-1,JJ-1)/DZ(JJ) 
TB(IJ,I1J2)=-RKLS*W0(J-1,JJ-1)/D2Z(JJ) 
ELSE 
TB(IJ,IJ1)=.5+RKLS*W0(J-1,J-1)/D2Z(JI) 
TB(IJ,IJ2)=.5—-RKLS*W0(J-1,J3J-1)/DZ(JJ) 
ENDIF 
CONTINUE 
DO 136 JJ=JIM+1,JIM+JE 
IJ3=(I-1) *JIM+JIJI 
TB(IJ,IJ3)=RHK*W0(J-1,JI3-1) 
CONTINUE 
TR(IJ)=RHS(J—-1) 
SWJ=0. 
DO 137 JJ=2,IM+JE 
SWJ=SWJ+W0(J-1,JJ-1) 
CONTINUE 
TR(IJ)=TR(IJ)+( TSLDO-TMPAIR) *(1.+RHK*SWJ ) 
ELSE 
TB(IJ,IJ)=1. 
TB(IJ,IJ-JM)=1. 
TR(IJ)=2.*(CT(J)-TMPAIR) 
ENDIF 
CONTINUE 


C....FOR PTS ON SOLID EXPOSED SURFACE 


155 


156 


DO 150 J=JM+1,IM+JE 

IJ=(I-1)*JM+J 

IF(NT.GT.NTCS) THEN 

TB(IJ,IJ)=1.+RHK*wW0(J-1,J-1) 

DO 155 JJ=2,9M 

IJ1l=(I-1) *IM+JIJI 

IJ2=(I-2) *IM+JJd 

TB(IJ,IJ1)=RKLS*W0(J-1,JJ-1)/D2Z(JI) 

TB(IJ,IJ2)#-RKLS*W0(J-1,J3J-1)/D2(Jd) 
CONTINUE 

DO 156 JJ=JM+1,JIM+JE 

IJ3=(I-1) *IM+JJI an 

IF(IJ3.NE.IJ)TB(IJ,IJ3)=RHK*W0(J-1,J3J-1) 
CONTINUE 
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TR(IJ)=RHS(J-1) 
SWwJ=0. 
DO 157 JJ=2,J5M+JE 
SWJ=SWJ+W0(J-1,JJ3-1) 
157 CONTINUE 
TR(IJ)=TR(IJ)+(TSLDO-TMPAIR) *(1.+RHK*SWJ ) 
ELSE 
TB(IJ,IJ)el. 
TR(IJ)=TSLDO-TMPAIR 


ENDIF 
150 CONTINUE 
C...FOR Jel. B.C. T(I,1)=T(1I,2) 
J=l 


DO 140 I#1,IM 
IJ=(I-1)*JM+d 
TB(IJ,IJ)el. 
TB(IJ,IJ+1)=-1. 


TR(IJ)=0. 
140 CONTINUE 
Cereic oetets tctate cbelela, cl chekdic ole cle lsielers else oe ese o.0.6.6 6 slelenel ofc sieleehaleie.e 6 66:6. 6 
CALL SOLVE(IM*JM+JE,TB,TR,TX,IWK,WK) 
Gaars 
KT=0 
DO 300 IJ#1,IM*JM+JE 
TN1(IJ)=TX(IJ) 
300 CONTINUE 
Gwe es 


CALL SURF(DTIME,IM,JM,RR,VO,AA,R,Z,DVDT, 
EWIS, COEFFH, SX, RHOW, CA) 
Craileretonetel o Sie as eG tohaleus ic Gceielernis eo anole sie baie eves tele tele leisihe ® ehete evahels © sie, 6 6 016 6 ee 
DO 417 KKeKNT,2,-1 
DO 417 IIel,JM+JE-1 
FTIM(II,KK)=FTIM(II,KK-1) 
417 CONTINUE 
DO 418 JJ=2,JM+JE 
J1l=(IM-1)*JM+JJ 
J2=(IM-2) *JM+JJI 
I=IM 
FTIM(JJ-1,1)=—-RHK*(TN1(J1)+TMPAIR-TSLDO ) 
418 CONTINUE 


i CALL NEXTR(NT,KNT, FTIM,RHS, TIME,W,F,WF,RSO) 

Ee IF (DVDT*DTIME.GT.VO)GO TO 3 

Costes eece rsa essaeeec acta eeeceeceecescceeececc seen ese eee eee ee: 

Cos eeee sec neeeeeeeeeec nec tecaeccaeec eee seca eeeeceeecereeeetens 
END 


C...SUB. ASSIGNING THE CENTER PTS. FOR THOSE CELLS AT EDGE WHICH ARE 
C...TRIANGLES. 
SUBROUTINE TRIO(IM,JM,R,Z,RO3,2Z03) 
DIMENSION R(IM+1,JM+1),Z2(IM+1,JM+1) 
DIMENSION RO3(IM) ,Z03(1IM) 
J=JM 
DO 10 I#1,IM 
A=(.5*(Z(1,J3)+2Z(I,J3+1) )-2(I+1,I9))/ 
$ (.5*(R(I,J)+R(1I,J+1) )-R(I+1,d) ) - 
B=a( .5*(Z2(I+1,5)4+Z(1I1,J+1) )-2(1I,75))/ 
$ (.5*(R(I41,9)4+R(1I,J3+1) )-R(I,Jd) ) 


83 


10 


RO3(I)=(Z(I,J)—-Z(I+1,J3)+A*R(I+1,J3)-B*R(I,J))/(A-B) 
Z03(I)=Z(1,J3)+B*(RO3(I)-R(I,J)) 

CONTINUE 

RETURN 

END 


Cure erence tateteretere 


SUBROUTINE NEWSO(KK,IM,JM,R,Z,RO3,Z203,1,J3,RN,2ZN,RE,ZE, 


Cc 


RW, 2ZW,RS,2ZS,RO,ZO) 


C....-SUBROUTINE FOR ASSIGNING THE POSITIONS AT NORTH, EAST, 
C...-CENTER OF A CELL. 
C....KK=4: CELL HAS FOUR SIDES; KK=3: CELL HAS THREE SIDES. 


C....THE LEFT UPPER CORNER IS THE POSITION R(I,J) AND 2(I,J) 


DIMENSION R(IM+1,JM+1),Z2(IM+1,JM+1) 
DIMENSION RO3(IM) ,Z03(IM) 
IF(J.EQ.JM) THEN 

RO=RO3 (I) 
20=2Z03(T) 
RN@RO3(I-1 
2ZN#=Z03(I-1 
RS#RO3(I+1 
Z2S8=203(I+1 
RWe.25*(R( 
ZWe .25*(Z( 
ELSE 
RO=.25*(R(I,J)+R(I1,J3+1)+R(I14+1,3+1)+R(I+1,J7) ) 
Z0=.25%*(Z2(1,9)4+2(1,53+1)4+2(I4+1,34+1)4Z2(I+1,7) ) 
RE=@=.25*(R(I,J+1)4+R(I,J3+2)4+R(141,3+2)+R(I+1,5+1) ) 
ZE=.25%*(2(1,5+1)42(1,5+2)42(I+1,3+2)+2(I+1,I5+1) ) 
RNe@=.25*(R(I,J+1)+R(I-1,J3+1)+R(I-1,J3)+R(I,J) ) 
ZN=.25*(Z(1,J3+1)+2Z(I-1,J3+1)4+2Z(I-1,J3)+2Z(1I,J) ) 
RW=.25*(R(I,J)4+R(I,J—-1)4+R(I4+1,J-1)+R(I+1,J) ) 
Z2W=.25*(Z2(1,3)4+2(1,3-1)+2Z(I+1,3-1)4+Z2(I+1,J9)) 
RS=.25*(R(I+1,J3)+R(I+1,J3+1)4+R(1+2,3+1)+R(I+2,d) ) 
ZS=@.25*(Z(I1+1,9)4+2(I4+1,5+1)4Z2(I4+2,5+1)4Z2(1I1+2,J9) ) 
ENDIF 

RETURN 

END 


) 
) 
) 
I,J)+#R(I,J-1)4R(I41,J-1)+R(I+1,3) ) 
I,J)4Z2(I,J—-1)4Z2(I4+1,J3—-1)4Z(I+1,79) ) 


SUBROUTINE RZ90(IM,JM,RR,AA,R,Z) 


»SUBROUTINE FOR ASSIGNING THE GRID POINTS..... 
-R(I,J), Z2(1,3): I=ml,IM+1, J=l,JM+1 
- IMAGINARY POINTS: 1) R(1,J), 2(1,d) 

2) R(IM+1,9), Z(IM+1,J) 

3) R(I,1), 2(1I,1) 


- AA: DROPLET HEIGHT; RR: DROPLET RADIUS. 


DIMENSION 2Z(IM+1,JM+1),R(IM+1,JM+1) 
RRR#=AA/RR 


DO 10 I=2,1M 
FIlel.-(1*1.-2.)/(IM*1.-2.) 
DO 10 J=2,JM 
GJ=(J*1.-2.)/(IM*1.-1.) 


IF(I.EQ.IM.OR.J.EQ.2) THEN 
IF(I.EQ.IM) THEN 

R(I,J)=#RR*GJ = 
Z(I,J)=0. 

ENDIF 


WEST, 


SOUTH 
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IF(J.EQ.2)THEN 
R(I,J)=0. 
Z(I1,/J)sAA* FI 
ENDIF 


ELSE 
FIR=1./(FI*RRR)-FI*RRR 
A=(1.+GJ**2+2.*(1.-GJ))*RR/FIR 
Be-2.*(GJ+(1.-GJ)/GJ)/FIR 
FACTOR=1. 
ELSE 
FACTOR#=1.E-6 
ENDIF 
ALFA1=(1.+B**2)*FACTOR 
BEITA=2.*(A*B-(GJ+(1.-GJ)/GJ)*RR) *FACTOR 
GAMMA=(A**2+(GJ**2+2.%*(1.-GJ))*RR**2) *FACTOR 
TF<BEITA**2-4. *ALFA1 *GAMMA 
Z(I1,J)=SQRT((1./(FI*RRR)+FI*RRR) **2/4.*RR**2-R(I,J) **2) 
$-(1./(FI*RRR)—-FI*RRR)/2.*RR 
10 CONTINUE 
DO 20 I#1,IM+1 
R(I,JM+1)=RR 
Z(I,JM+1)=#0. 
20 CONTINUE 
C..FOR IMAGINARY POINTS 
DO 30 I=#1,IM+1 
R(I,1)=#-R(I,3) 
$+2.*R(I,2) 
Z(I,1)=Z2(1I,3) 
30 CONTINUE 
DO 40 J=1,JM 
R(IM+1,J)=R(IM-1,J) 
Z(IM+1,J)=-Z(IM-1,J) 
$+2.*Z(IM,J) 
40 CONTINUE 
I=1 
AA1L@=AA+(Z(2,2)-2(3,2) ) 
RRR=AA1/RR 
FIel. 
Z(I,2)=FI*AA1L 
DO 50 J=3,JM 
GJ=(J*1.-2.)/(JIM*1.-1.) 
FIR=1./(FI*RRR)-FI*RRR 
A=(1.4+GJ**2+2.*(1.-GJ))*RR/FIR 
B=-2.*(GJ+(1.-GJ)/GJ)/FIR 
IF(FIR.GT.1.E-3) THEN 
FACTOR=1. 
ELSE 
FACTOR=1 .E-6 
ENDIF 
ALFA1=(1.+B**2)*FACTOR 
BEITA=2.*(A*B-(GJ+(1.-GJ)/GJ)*RR) *FACTOR 
GAMMA=(A**2+(GJ**2+2.*(1.-GJ))*RR**2) *FACTOR 
TF=BEITA**2—4. *ALFA1*GAMMA 
IF(TF.LT.0.)TF=#0. 
Z(I,J)=SQRT((1./(FI*RRR)+FI*RRR) **2/4.*RR**2-R(I,J) **2) 
$-(1./(FI*RRR)-FI*RRR)/2.*RR 
50 CONTINUE 
R(1,1)"#-R(1,3)4+2.*R(1,2) 
Z(1,1)=2(1,3) 


Crarenens 


RETURN 

END 
Cc 

SUBROUTINE DZ0(1IM,JM,R,Z,RO3,203,DZ) 
Cc 


C...-eSUBROUTINE FOR FINDING THE DZ 
DIMENSION R(IM+1,IM+1),Z(IM+1,JIM+1) ,DZ(JM) 
DIMENSION RO3(IM),2Z03(IM) 
I=@=IM 
DO 135 JJ=2,JM 
IF(JJ.EQ.JM) THEN 
RM=RO3(TI) 
ZM=203(1) 
RP=RO3(I-1) 
ZP=Z03(I-1) 
ELSE 
RM=.25*(R(I,JJ)+R(I,JI+1)+R(I+1,J3J+1)+R(I+1,J5J) ) 
ZM@=.25%*(Z2(1,33)+Z(1,JJ+1)4+Z2(I+1,JJ+1)+Z(I+1,5I) ) 
RP=.25*(R(I-1,JJ)+R(I-1,IJ+1)4R(I,JJ+1)4+R(I,J3J) ) 
ZP=.25*(Z(I-1,5I)42(I-1,JI+1)4+2(1,35+1)4+Z2(1I,5I) ) 


ENDIF 
DZ(JJ)=SQRT( (RM-RP) **2+(ZM-ZP) **2) 
135 CONTINUE 
RETURN 
END 
Cc 
SUBROUTINE SOLVE(N,A,B,X,IWK,WK) 
Cc 


C....-SUBROUTINE FOR SOLVING MATRIX 

C...’LINK XXX,SOLVE,IMSL/LIB’.... 
DIMENSION A(N,N),B(N),X(N) 
DIMENSION IWK(N) ,WK(N+N*(N-1)/2) 
INTEGER N 3 
CALL L2NRG(N,A,N,A,N,WK, IWK) 
DO 100 J=#1,N 
X(J)=0.0 
DO 100 I=#1,N 
X(J)=X(J)+A(J,1)*B(I) 

100 CONTINUE 


RETURN 
END 
c 
SUBROUTINE SURF(DTIME,IM,JM,RR,VO,AA,R,Z,DVDT, 
$ EWIS, COEFFH,SX,RHOW,CA) 
Cc 


C....-SUBROUTINE FOR CALCULATING THE SIZE OF THE LIQUID DROP RR, AA 
C....-WHILE RR IS A CONSTANT 
Cc 


DIMENSION R(1+IM,1+J3M) ,2(1+IM,1+JM) ,EWIS(JM),SX(JM) 
COMMON /MATER/ LMAT 
C...REMEMBER TN1+TEMPO IS THE REAL TEMP. FOR EVALUATING THE PROPERTIES. 
PI=4.*ATAN(1.) 
SSX=0. 
DO 10 J=2,JM 
RJ=.5*(R(2,5)+R(2,J5+1) ) 
DRJ=SQRT((R(2,0+1)—-R(2,0) ) **2+(2(2,5+1)-2(2,59) )**2) 
SSX=SSX+EWIS(J) *COEFFH*SX(J)*RJ*DRJ © 
10 CONTINUE - 
DVDT#=2.*PI*.624/( RHOW*CA) *SSX 
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Cc 


VO=VO-DVDT*DTIME 

BETA#=RR/EXP(1./3.*LOG(.75*VO/PI) ) 

Pl=1./2. 

P2=1./3. 
UsEXP(1./3.*LOG(3.*VO/PI+SQRT(RR**6+9./PI**2*VO**?) ) ) 
AA®U-RR**2/U 

RETURN 

END 


149 


SUBROUTINE NEXTR(NT,KT,FTIM,R,TIME,W,F,WF,RSO) 


DIMENSION W(NP,NP),F(NP),FTIM(NP,KT) 
DIMENSION WF(NP),R(NP),RSO(NP+1) 
COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NDP,NP,TN 
COMMON /INTG/ ATO,ERRREL, IRULE 
COMMON /MATER/ LMAT 
DO 149 II=1,NP 
R(II)=0. 
CONTINUE 
IF(NT.LE.KT) THEN 
KKT=NT 
ELSE 
KKT=KT 
ENDIF 
DO 150 IT=1,KKT 


SELECT PROPER WEIGHTS 


C.. 


Cc 


Cc 
Cc 


PICK 


100 


TO=IT*DTO 


CALL WEIGHT(TO,W,RSO) 
UP CORRESPONDING FLUX 
DO 100 I=1,NP 
F(I)=FTIM(I,IT) 
CONTINUE 


SUMMATION IN SPACE 


CALL MULT(NP,NP,W,F,WF) 


SUMMATION IN TIME 


150 


50 


CALL ADD(NP,WF,R) 
CONTINUE 

RETURN 

END 


SUBROUTINE MULT(M,N,W,F,WF) 
REAL W(M,N),F(N) ,WF(N) 
DO 50 I=1,M 

WF(I)=0.0 . 

DO 50 K=1,N 
WFE(I)=WF(I)+W(I,K)*F(K) 
CONTINUE 

RETURN 

END 

SUBROUTINE ADD(NP,WF,R) 
REAL WF(NP) ,R(NP) 


DO 50 I=#1,NP 
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R(I)=#R(1I)+WF(I) 
50 CONTINUE 

RETURN 

END 


SUBROUTINE WEIGHT(T0O,W,RSO) 


REAL W(NP,NP),RSO(NP+1) 
REAL QDAG,BSIOE,ERF 
EXTERNAL Fl 


COMMON /CONST/ ALFA, PI 

COMMON /GRID/ DRO,DT0,NDP,NP,TN 
07800 

COMMON /INTG/ ATO,ERRREL, IRULE 

COMMON /RRO/ R,RO 

COMMON /MATER/ LMAT 


DO 50 Ie#1l,NP 

DO 50 J=#1,NP 
DRO=-RSO(J+1)-RS0(J) 
R=.5*(RSO(I+1)+RS0(I 
RO=.5*(RSO(J+1)+RS0( 
® SET UP (W) FOR TO < DTO 
ERRABS=0.0 
IF(TO.LT.DTO) THEN 
W(I,J)=QDAG(F1,AT0,0.5,ERRABS, ERRREL, IRULE, RESULT, ERREST) 
W(I,J)=RESULT 

END IF 


)) 
J)) 


» 


SET UP (W) FOR TO > DTO 
IF(TO.GE.DTO) THEN 
ARG1=DR0/SQRT(16*ALFA*TO ) 
ARG3=(R-RO) **2/( 4*ALFA*TO) 
ARG4=0.5*R*R0/(ALFA*TO ) 
IF(I.NE.J) THEN 
W(I,J)=(1.0/SQRT(4.0*PI*ALFA) )*RO*(TO**(-1.5))* 
& BSIOE(ARG4) *EXP(-ARG3 ) *DRO*DTO 
ELSE 
W(I,J)=(R/TO) *BSIOE(ARG4) *ERF(ARG1) *DTO 
END IF 
END IF 


50 CONTINUE 
RETURN 
END 


FUNCTION F1(TO) 

REAL BSIOE,ERF 

COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NP 
COMMON /RRO/ R,RO - 


ARG1=DRO/SQRT(16*ALFA*TO ) 
ARG3=(R-RO) **2/( 4*ALFA*TO ) 

ag Sie beh oe 

E=0.0 

IF(ARG3.LE.80.0) E=EXP(-ARG3) 

IF(R.NE.RO) THEN 

Fl=(1.0/SQRT(4.0*PI*ALFA) )*RO*(TO**(-1.5))* 


& BSIOE(ARG4) *E*DRO 
ELSE 
Fl=(R/TO) *BSIOE(ARG4) *ERF(ARG1 ) 
END IF 
RETURN 
END 
Cc 
SUBROUTINE CTEMPO(IM,JM,JE,TREF,TWATER, TSOLID,CT) 
Cc 
C...SUBROUTINE CALCULATING THE CONTACT TEMPERATURE TO BE USED AT 
C...THE SOLID-LIQUID INTERFACE 
DIMENSION TWATER(IM*JM+JE) ,CT(JM) 
COMMON /SOLID/ SK,SC,SR 
COMMON /WATER/ COEFKW, CPW, RHOW 


Co cecce 
SQRCKW=SQRT( RHOW*CPW* COEFRW ) 
SQRCKA=SQRT(SR*SC*SK) 
Cieccets 
DO 10 J=1,JM 
CT(J)=(SQRCKW* ( TWATER( (IM—2) *JM+J)+TREF ) 
$ +SQRCKA* TSOLID) /( SQRCKW+SQRCKA ) 
10 CONTINUE 
RETURN 
END 


C....-SUB. FOR ASSIGNING POINTS ON SOLID SURFACE 
SUBROUTINE RSD(JM,JS,RR,RSO) : 
DIMENSION RSO(JM+4*JS) 
DR1=RR/(1.*JM-—-1. ) 

DO 10 J=1,JIM 
RSO(J)=(J-1.)*DR1 

10 CONTINUE 
RSO(JM+1)=#@RSO(JM)+DR1 
DO 20 J=JM+2,IM+JS 
RSO(J)=RSO(J-1)+DR1 

20 CONTINUE 
DO 30 J=JIM+JS+1,IM+2*IS 
RSO(J)=RSO(J-1)+2.*DR1 

30 CONTINUE 
DO 40 J=JM+2*IS+1,IM+3*JS 
RSO(J)=RSO(J-1)+4.*DR1 

40 CONTINUE 
DO 50 J=JM+3*JS+1,IM+4*JIS 
RSO(J)=RSO(J-1)+8.*DR1 

50 CONTINUE 
RETURN 
END 
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